Relaxation dynamics of the Holstein polaron 
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Keeping the full quantum nature of the problem we compute the relaxation time of the Holstein 
polaron after it was driven far from the equilibrium by a strong oscillatory pulse. Just after the 
pulse the polaron's kinetic energy increases and subsequently exhibits relaxation type decrease with 
simultaneous emission of phonons. In the weak coupling regime partial tunneling of the electron from 
the polaron self-potential is observed. The inverse relaxation time is for small values of electron- 
phonon coupling A linear with A, while it deviates downwards from the linear regime at A > 0.1/wo- 
The imaginary part of the equilibrium self energy shows good agreement with the inverse relaxation 
time obtained from nonequilibrium simulations. 

PACS numbers: 71.38.-k, 63.20.kd, 72.10.Ht 
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Photocxcitation is due to recent technological advances 
in ultrafast spectroscopy becoming one of the main ex- 
perimental approaches to disentangle different elemen- 
tary excitations in a real-time domain on a femtosecond 
scale. This rapidly developing field enables a novel in- 
sight into quantum many-body systems, however it also 
requires the development of new theoretical concepts. 
In this context, numerical simulations of nonequilibrium 
quantum mechanical systems may provide a key insight 
into processes on a femtosecond time scale. 

Several pump-probe experiments reported the observa- 
tion of self-trapping of excitions emerging after the pump 
pulse [lj-Q • These experiments then triggered theoretical 
studies of different scenarios of the polaron formation, re- 
vealing a complex interplay between a single electron and 
quantum phonons under nonequilibrium conditions 0, [B| . 

The role of electron-phonon (e-ph) interaction in sev- 
eral different classes of strongly correlated materials is 
despite intensive research still ambiguous since a subtle 
interplay between electron-electron and electron-phonon 
interaction may lead to various unconventional proper- 
ties. Pump-probe techniques have shown a potential to 
identify fingerprints of these interactions during the re- 
laxation process ■ One of the main challenges of such 
experiments represents the choice of underlying theoret- 
ical framework to interpret results. Recently, the phe- 
nomenological two temperature model developed long 
time ago jlClL EH , was extended to more involved ap- 
proaches like the three-temperature model and Boltz- 
mann equation approaches (l3l - [l5j . Moreover, the as- 
sumption of fast relaxation within the electron subsystem 
was recently challenged in the case of strong electronic 
correlations (l6| . In contrast to many previous findings 
it was shown that in the ID strongly correlated system 
coupled to phonons, the relaxation on short time scales 
is mostly due to e-ph interaction [l7j . 

Theoretical studies of the polaron motion in strong but 



constant electric field started with the seminal work by 
Thornbcr and Feynman fl8j . Later works mostly relied 
on the rate or Boltzmann equations U 21 \ ■ While quan- 
tum coherent effects arc absent in the Boltzmann de- 



scription |22j , some recent approaches [23|, |2J] show that 
taking fully into account quantum effects is decisive to 
obtain proper electric field dependence of the drift ve- 
locity at large electric fields. Since most past works fo- 
cused on the influence of a constant electric field on the 
Holstein polaron, the impact of a short oscillatory pulse 
on polaron relaxation dynamics remains an open prob- 
lem despite significant advances in ultrafast spectroscopy 

BEES. 

In this Letter we present results of a fully quantum 
mechanical time evolution of the Holstein model driven 
far from the equilibrium by a laser pulse. We determine 
characteristic relaxation times for the system which is 
initially in equilibrium at zero temperature. The investi- 
gated Holstein polaron is subjected to a spatially homo- 
geneous and time-dependent scalar potential that mimics 
a short laser pulse: 



<f>(t) = Ae-W-^M* wn(u p (t-t c )), 



(1) 



which is incorporated into the Hamiltonian via a Peierls 
substitution in the hopping amplitude: 



H = -t, 



(2) 



where cj and a\ represent electron and phonon creation 

operators at site i, and rii = c|cj is the electron density, 
wo denotes the dispersionless phonon frequency and to is 
the nearest neighbor hopping amplitude. 

The system is described by two dimensionless param- 

g 2 /(2u t ) that 



eters a and A where a 



¥>■ and A 

to 



2 



determine the crossover from adiabatic(a < 1) to nona- 
diabatic (a ^> 1) limit and the weakly dressed electron 
(A < 1) to a heavy polaron (A 3> 1)- We measure the 
electric field F = —d t 4>{t) in units of to/eoa, where eo is 
the unit charge and a is the lattice distance. The unit of 
energy is given by hopping amplitude to and the unit of 
time is h/to- From here on we set a = eo = H = to = 1. 
Unless explicitly stated the phonon frequency is set to 
wo = 1- The pulse in Eq. (JTJ is centered at t c = 5, while 
the width is given by td = 2. 

We solve the time-dependent Schrodingcr equation for 
a single Holstein polaron on an infinite one-dimensional 
chain. We use the numerical method based on the ex- 
act diagonalization of the variational Hilbcrt space that 
led to numerically exact solutions of the polaron ground 



27 1 and low- lying excited-state properties 28l 



state 

as well as for description of the time-dependent case 
The total energy gain from the external pulse is given by 
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A(H(t)) = \ (j(t))F(t)dt, 



(3) 



where j(t) = —d<j,H is the current operator. Equation (j3]) 
is as well used as an additional time propagation accuracy 
check, see also Ref . (30| . 

Before examining the more physically relevant case we 
analyze the action of a pulse in a form of a delta function 
F(t) = —$o5(t) that leads to a simple form of the scalar 
potential <j>(t) = <&o9(t) where 0(t) is the Heaviside func- 
tion. Action with 4>(t) on a free electron state (A = 0) at 
k = shifts its kinetic energy from E^ = — 2 at t < 
to .Ekm = — 2cos($o) when t > 0, leaving the electron 
in an excited, but an eigenstate at t > 0. In Fig. QJa) 
we show results for A = 0.1 and $o = i"- Just after the 
pulse, the increase of the total AE(t) = (H(t)) — (H(t = 
0)), as well as the kinetic energy AE^i n (t), reach value 
AE ~ AE kin — -2cos($ ) +2 = 4 while the change of 
the phonon energy AE pn remains close to its value in the 
polaron ground state, i.e. AE p h ~ 0. After initial time 
t ^ t% ~ 10, AE\<_i n (t) exhibits a relaxation type expo- 
nential decay towards a constant value AE\an(t — ► oo)j 
meanwhile AE p h(t) increases and the electron-phonon 
interaction term Ai? e _ p h(t) remains nearly a constant. 
This dynamics is interpreted as a transfer of the ex- 
cited electron kinetic energy into phonon excitations. In 
Fig. QJb) only AE\^ n (i) is shown for different choices of 
pulse amplitude $o- In ah cases except for $o = 7r/4, 
AEk in (t) decrease exponentially with roughly the same 
relaxation time. In the case of $o = 7r /4 the increase 
of AEki n (t) is lower than ujq indicating that the energy 
transfer from the pulse was insufficient to allow electron 
relaxation via phonon emission, thus no relaxation is ob- 
served. 

We next consider a more realistic form of a pulse de- 
scribed by the scalar potential of Eq. ([TJ. At small 
A = 0.1 we observe a gain of the total energy AE(i), see 
Fig. [He), signaling that despite rather weak A, the sys- 



tem has absorbed a substantial amount of energy. Note 
that after the pulse is switched off for t > 9 = t D g, 
the total energy remains constant while there is a clear 
redistribution between expectation values of the parts 
of the Hamiltonian. Redistribution between the kinetic 
AE\ti n (t) and the phonon AE pn (£) part of the total en- 
ergy clearly indicates a relaxation of the system. After 
the pulse, AE p h(t) increases with time and AE^m (t) ex- 
ponentially decreases, while AE e - p h(t) oscillates roughly 
around zero. This behavior is similar to the case when 
4>(t) = ®o@(t), see Fig. [JJa). For t > t g most of the 
gained energy is absorbed by the lattice. We should re- 
mark that the relaxed expectation values of the kinetic 
energy are not the same as before the pulse. As a gen- 
eral rule we find in all other cases AE]^ n (t — > oo) < ujq, 
which indicates that a finite value of AE k i n (t — > oo) is 
a consequence of the gap for optical phonons [22j]. Nev- 
ertheless, we detect a clear tendency of the exponential 
decay of the kinetic energy towards only a slightly ele- 
vated value AEkin (i: — > oo) in comparison to the initial 
energy E kin (t = 0). We have also computed the relax- 
ation of the kinetic energy for different amplitudes of the 
pulse, presented in Fig. [TJd). We found that the decay 
time is within our numerical accuracy independent of the 
amplitude of the pulse A as long as A exceeds a threshold 
value A > 1, see also the inset of Fig. |Tf d) , which is in 
agreement with the Boltzmann theory. 

Our results are qualitatively consistent with experi- 
ments on quasi-one-dimensional systems j3lj . where ex- 
ponential relaxation takes place [l|, [2| . Relaxation was 
also observed in 2D systems 26|, |32|, |33(, where it is a 
consequence of the molecular dipoles rearrangement in 
the vicinity of the electron, which can be described by 
formally equivalent electronic polaron model (34| . 

In the intermediate coupling regime, namely A = 1.0, 
a different response of the system is expected due to a 
bigger gap between the polaron band and the continuum 
of excited states. After the pulse almost the entire ex- 
cess energy is absorbed into the lattice vibrations, i.e. 
AE(t > t oS ) — AE ph (t > t oS ), see Fig.QJe). In addition 
to a large increase, the latter displays as well oscillations 
with the period corresponding roughly to the phonon fre- 
quency. In Fig.[]Tf) we present results of A£^in(*) f° r dif- 
ferent amplitudes of the pulse. Apart from pronounced 
oscillations, we again observe approximately exponential 
decay in the kinetic energy, roughly independent of the 
strength of A when A > 1. We also observe distinct 
long-time limits of kinetic energies, however in all cases 

AEkin (t > toff) < wo- 

To get further insight into the relaxation dynamics we 
calculated the average number of phonon quanta located 
at a given distance x from the electron 

70) = (^rna\ +x a i+x ), (4) 

i 

satisfying the sum rule (n p h) = ^2 x l( x )- At time t = 
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FIG. 1: (a) Expectation values of different parts of the Hamil- 
tonian vs time at A = 0.1 and uo = 1.0, after action of 
4>(t) = $ o 0(t), for $ = 7T. (b) A£ kin (i) vs t for different 
$o = 7r, 37r/4, 7r/2, 7r/4, while other parameters are the same 
as in (a), (c) Expectation values of different parts of the 
Hamiltonian vs time at A = 0.1 and ljo = 1, after action of 
pulse as given in Eq. {TJ with A = 3, where ui p — 1.5 is set 
close to the maximum of the optical conductivity, while pa- 
rameters t c — 5 and t p — 2 remain unchanged throughout this 
work, (d) ASkin(t) vs t for different values of A = 1, 3, 5, the 
rest is the same as in (c). The inset represent the difference 
between kinetic energy and E^ in {t — ► oo) for the same pulse 
amplitudes as the main figure, (e) Energy of the different 
parts of the Hamiltonian vs t for A = 1.0 and uj p = 2.5. (f) 
A£/ki n (t) vs t for different values of A = 1, 3, 5. while the rest 
is the same as in (e). 



0, 7(x) displays a pronounced peak at the position of 
the electron (x = 0), consistent with the shape of the 
polaron in equilibrium. In the weak coupling regime, at 
A = 0.1 and t > t s, j(x) shows beside original polaron 
correlation peak also pronounced peaks separating from 
the central peak in both directions as time increases, see 
Fig. HJa). This result is consistent with a hypothesis 
that a strong pulse splits the polaron into an excited 
polaron and a nearly free electron. The wavefunction is 
a superposition of an excited polaron, responsible for the 
large value of j(x) at x — and a nearly free electron, 
traveling predominantly in the x > direction. The 
more pronounced signal for 7(2; < 0) can be interpreted 
as a partial tunneling of electron part of the wavefunction 
from polaron self-potential that remains located at x — 0. 
The asymmetry on the parity transformation x — > — x 
is dynamically induced and can be tuned by choosing 
different shape of the incoming pulse. 

The peak at x < starts to diminish with time be- 
cause: (i) the escaped nearly free electron is gradually 
captured by the lattice and (ii) the excited electron redis- 
tributes its excess energy into constantly spreading area 
of excited lattice vibrations giving rise to nearly uniform 



FIG. 2: Electron-phonon correlation function 7(1) for pulse 
amplitude A = 3; (a) A = 0.1 and the pulse frequency u) p = 
1.5 and (b) A = 1.0 and uj p — 2.5, at different times. Note 
that there are different scales on both figures. 



but elevated values of j(x), clearly seen in Fig. HJa). 
Here we should stress that redistribution in the corre- 
lation function must be due to electron's motion since 
phonons are dispersionless. This hypothesis is well sup- 
ported by the estimation of the velocity of the side peaks 
v p = Ax/ At ~ 2 representing the maximal group ve- 
locity of the weakly coupled electron [see dashed line in 
Fig. [2Ja)] . Similar partial tunneling was noticed within 
the adiabatic limit of the driven Su-Schricffcr-Hccgcr 
problem (35j . Above considerations give complemen- 
tary real space interpretation of related experimental re- 
sults, where excitations are identified in the frequency 
domain [l|, In the strong coupling regime, namely 
A = 1.0 as presented in Fig. [2jb), the polaron peak at 
7(2; ~ 0) is preserved, but it broadens with time and we 
observe no peak due to partial electron tunneling. 

We computed the relaxation time by fitting the expec- 
tation value of the kinetic energy after the pulse with the 
simple expression A£^ in (t) = AE^ hl (t — s- 00) + Be~ l l T , 
where A£^i n (i — s- 00) is the kinetic energy after relax- 
ation and t is the relaxation time, see Fig. [31 The 
inverse relaxation time 1/r in the extreme weak cou- 
pling A < 0.1 regime shows linear increase with electron- 
phonon coupling, consistent with 1/r w 2loq\. Consid- 
ering an emission of a phonon by the excited electron 
with the kinetic energy above onc-phonon threshold us- 
ing Fermi golden rule yields linear dependence on A, i.e. 
1/r = 2woA/ sin(fc/), where kf is the final electron's mo- 
mentum and the inverse relaxation time is determined 
by the longest decay time, namely at kf = 7r/2. With 
increasing A the inverse relaxation time 1/r first devi- 
ates downwards from the linear A dependence and then 
saturates as it becomes comparable to the pulse width. 



4 



0.35 
0.3 
0.25 

: 0.2 

0.15 
0.1 
0.05 



- Gauss. 
2A(d 
-lm[Z] 

* e(t). 



c =1.0 




0.1 0.2 0.3 

A (d n 



0.1 



0.2 0.3 

A 



0.4 



0.5 



FIG. 3: Inverse relaxation time 1 /r (squares) vs A for Uo =1.0 
using Gaussian pulse, Eq. ([T]). The dashed line represents 
1/r = 2woA, see discussion in text. Diamonds: the imagi- 
nary part of the averaged equilibrium self energy Im[E]. Cir- 
cles: 1/r after the instantaneous pulse <f>(t) = $o6(t) (circles) 
where pulse strength $o = w was used. The horizontal dot- 
ted line indicates the inverse pulse width. The inset presents 
1/r vs. Xloq, using Gaussian pulse, Eq. JTJ), f° r wo = 1 (full 
squares) and luq —0.7 (open squares). 



The larger error bars are a consequence of strong oscilla- 
tions and a smaller spatial extend of the variational phase 
space in the strong coupling regime. Calculating 1/r at 
smaller cjo = 0.7 shows that the scaling of 1/r ~ /(Awrj), 
where f(x) is some unknown function [note that within 
Fermi golden rule f(x) = 2x], persists beyond the linear 
in A regime, see the inset of Fig. [3] The saturation of 1/r 
is absent in the case when <fi(t) = §>o9(t), also presented 
in Fig. [3] The relaxation time after an instantaneous 
pulse corresponds to the process averaged over all fre- 
quencies and this explains the deviation of the inverse 
relaxation time from a finite width pulse. 

An alternative method for the computation of 1/r 
from equilibrium properties is via the imaginary part 
of the self energy Im[E(u))] [36| that represents the in- 
verse of the relaxation time of the quasiparticle excitation 
1/rs = Im[£(tjp)]. Since Im[S(w)] depends on the value 
of the frequency of the pulse u p we define 1 /ts via the av- 
erage value t(w p ) = (1/V27Tct) / C^-^/W^^da}, 
where a = 1/td, and integrate over a distribution of fre- 
quencies corresponding to the Fourier transform of the 
pulse. Although there is no a priori reason that this equi- 
librium quantity represents the inverse relaxation time 
even after the strong non-equilibrium process the agree- 
ment with the actual values of the inverse time is excel- 
lent in the weak coupling limit. 

To conclude, in this Letter we studied the relaxation 
dynamics of the Holstein polaron after the strong photo- 
excitation. In all cases a threshold value of the absorbed 
energy AE ~ ujq exists above which the relaxation dy- 
namics via phonon emission is observed. We computed 
the relaxation time that is mostly independent of the 
shape of the pulse. We focused on the values of e-ph 



coupling A below the crossover to the strong-coupling 
regime, i.e. A c ~ 1. In this range of A, relaxation dy- 
namics exhibits two distinct regimes with qualitatively 
different behavior: the regime of very weak c-ph cou- 
pling, A < 0.1/uo, and the regime when A > 0.1/ujq. 

In the weak-coupling regime A < 0.1/wo, 1/r roughly 
follows the linear scaling 1/r ~ 2woA obtained from the 
Fermi golden rule. The relaxation process is described 
by the lowering of the kinetic energy at the expense of 
the increased phonon energy. A partial tunneling of the 
electron from the polaron self-potential occurs, result- 
ing in a rapid increase in time of the average distance 
between the electron and lattice deformation, followed 
by the subsequent re-trapping process. Moreover, values 
of r obtained from real-time calculations can be as well 
reproduced by calculation of the imaginary part of the 
averaged electron self-energy. 

In the regime of A > 0.1 /ujq, 1/r deviates from the 
linear in A dependence well below the crossover to small 
polaron regime. Most of the energy absorbed by the sys- 
tem is again deposited into phonon excitations, however, 
in contrast to A < 1, phonon excitations remain in the 
close proximity of the polaron. The real-time calculation 
reveals oscillations in Ai?km (t) and other expectation val- 
ues, with the period of the phonon frequency, T ~ 2tt/u}q. 
This result is in agreement with a recent study of a half- 
filled 2D Hubbard-Holstein model [33] , where oscillations 
with the period T were observed for A < A c . 

In comparison with the experiments our relaxation 
times are short, which is a consequence of rather high adi- 
abatic coefficient a = 1. While in the regime A < 0.1/wo 
the effect on r by lowering a can be obtained from the 
Fermi golden rule r oc 1/Awrj = 1/ Actio, our numerical 
results show that even in the regime when A > 0.1 /ujq 
there exists approximate scaling r ~ l//(Awo) that can 
be used to extend our results towards potentially more 
physically relevant values of a and consequently longer 
relaxation times r. Another interesting topic for further 
research is the effect of the dimensionality, since our sim- 
ulations agree with the exponential relaxation in quasi 
ID system, while relaxation is slower in 2D experiments. 

Let us briefly discuss the possible relevance of our re- 
sults to correlated electron systems with a finite elec- 
tron density, such as e.g. the Hubbard model. Re- 
cent noncquilibrium DMFT study showed that the re- 
laxation time of the pump-excited Hubbard model in the 
case of large Coulomb repulsion is unexepctedly long [l6| 
due to exponentially slow decay rate of pump-generated 
doublons, in agreement with experiments on optical lat- 
tices HU. These results open a relevant question about 
the dominant mechanism of fast relaxation observed in 
photoexcited strongly correlated materials. Lately, the 
emission of phonons in the ID Hubbard-Holstein model 
was indeed shown to be a very efficient relaxation mech- 
anism 17J , where a non negligible amount of phonons is 



already emitted during the application of the pulse. This 



5 



observation is in agreement with our results presented in 
Fig. [I] and suggests that the time-evolution within the 
full nonequilibrium process [both during and after the 
pulse] may provide a comprehensive understanding of the 
photoexcited polaronic systems. 
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